Chapter 4 Diversity analysis

load("data/data.Rdata")
load("data/gifts.Rdata")
load("data/beta_19022025.Rdata")

genome_counts_filt <- genome_counts_filt %>% 
  column_to_rownames(var = "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0)%>% 
  select(where(~ sum(.) > 0)) %>% 
  select(-AC85) %>% 
  rownames_to_column(., "genome")

genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)

table <- genome_counts_filt %>%
  t() %>%
  row_to_names(row_number = 1) %>% 
  as.data.frame() %>% 
  mutate_if(is.character,as.numeric) %>% 
  rownames_to_column(., "sample")
  
master_table <- sample_metadata %>%
  mutate(sample=Tube_code) %>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a")) %>% 
  filter(Tube_code %in% table$sample) %>%
  mutate(treatment = sub("^\\d+_", "", time_point))%>% 
  left_join(., table, by=join_by("Tube_code" =="sample"))

sample_metadata <- master_table %>% 
  select(1:13)

genome_counts_filt <- master_table %>% 
  select(12,14:323) %>% 
  t() %>%
  row_to_names(row_number = 1) %>% 
  as.data.frame() %>% 
  mutate_if(is.character,as.numeric) %>% 
  filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
  dplyr::select(where(~ !all(. == 0)))%>% 
  rownames_to_column(., "genome")

genome_counts_filtering <- master_table %>% 
  select(12,14:323) %>% 
  t() %>%
  row_to_names(row_number = 1) %>% 
  as.data.frame() %>% 
  mutate_if(is.character,as.numeric) %>% 
  filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
  dplyr::select(where(~ !all(. == 0)))

genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
match_data(genome_counts_filtering,genome_tree)

genome_metadata <- genome_metadata %>% 
  filter(genome %in% genome_counts_filt$genome)

4.1 Calculate diversities

4.1.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]

genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]
genome_counts_filt1 <- genome_counts_filt1 %>% 
  remove_rownames() %>% 
  column_to_rownames(var = "genome") %>% 
      filter(rowSums(. != 0, na.rm = TRUE) > 0) %>% 
  rownames_to_column(., "genome")

dist <- genome_gifts1 %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt1 %>%
  filter(genome %in% rownames(dist)) %>% 
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample))%>%
  full_join(functional, by = join_by(sample == sample))

treatment_colors <- c('#008080', "#d57d2c")

4.1.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]

dist <- genome_gifts1 %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

beta_q1f <- genome_counts_filt1 %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

save(beta_q0n, 
     beta_q1n, 
     beta_q1p, 
     beta_q1f, 
     file = "data/beta_19022025.Rdata")

4.2 Does the antibiotic administration change the microbial community?

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") ) 

4.2.1 Alpha diversity

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") ) %>% 
  ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(neutral ~ time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(4.0048)  ( log )
Formula: neutral ~ time_point + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   223.7    229.6   -107.9    215.7       28 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4630 -0.6452 -0.1754  0.6566  2.6564 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.07822  0.2797  
Number of obs: 32, groups:  individual, 18

Fixed effects:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)               2.9152     0.1558   18.72  < 2e-16 ***
time_point2_Antibiotics  -0.9432     0.2144   -4.40 1.08e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
tm_pnt2_Ant -0.538

Neutral

Model_neutral <- lme(fixed = neutral ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
      AIC      BIC    logLik
  225.017 230.6218 -108.5085

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.238645 7.593502

Fixed effects:  neutral ~ time_point 
                            Value Std.Error DF   t-value p-value
(Intercept)              19.29561  2.000903 17  9.643449   0e+00
time_point2_Antibiotics -11.58649  2.716098 13 -4.265859   9e-04
 Correlation: 
                        (Intr)
time_point2_Antibiotics -0.631

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.86139926 -0.51381109 -0.09492589  0.56226694  1.96603543 

Number of Observations: 32
Number of Groups: 18 

Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  138.9014 144.5062 -65.45072

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    1.046827 1.691404

Fixed effects:  phylogenetic ~ time_point 
                            Value Std.Error DF   t-value p-value
(Intercept)              5.461590 0.4814205 17 11.344740  0.0000
time_point2_Antibiotics -0.811302 0.6097318 13 -1.330588  0.2062
 Correlation: 
                        (Intr)
time_point2_Antibiotics -0.584

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.65239178 -0.50119332 -0.03450314  0.42332882  1.59503693 

Number of Observations: 32
Number of Groups: 18 

Functional

Model_funct <- lme(fixed = functional ~ time_point, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_funct)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
        AIC       BIC   logLik
  -14.95179 -9.346998 11.47589

Random effects:
 Formula: ~1 | individual
         (Intercept)  Residual
StdDev: 1.203543e-06 0.1504954

Fixed effects:  functional ~ time_point 
                             Value  Std.Error DF  t-value p-value
(Intercept)              1.4855915 0.03650050 17 40.70058  0.0000
time_point2_Antibiotics -0.0344778 0.05331239 13 -0.64671  0.5291
 Correlation: 
                        (Intr)
time_point2_Antibiotics -0.685

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.2430042 -0.5618554  0.2917987  0.6160260  2.0402200 

Number of Observations: 32
Number of Groups: 18 

4.2.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics")) %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.013523 0.0135230 2.0312    999   0.17
Residuals 30 0.199726 0.0066575                     

adonis2(richness ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.575038	0.1307754	4.513519	0.001
Residual	30	10.468804	0.8692246	NA	NA
Total	31	12.043842	1.0000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03475 0.034752 2.5153    999  0.105
Residuals 30 0.41450 0.013817                     

adonis2(neutral ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.778963	0.1605587	5.738054	0.001
Residual	30	9.300869	0.8394413	NA	NA
Total	31	11.079832	1.0000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.31299 0.312989 16.486    999  0.002 **
Residuals 30 0.56954 0.018985                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(phylo ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.569966	0.2885492	12.16735	0.001
Residual	30	3.870929	0.7114508	NA	NA
Total	31	5.440895	1.0000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.18537 0.185373 4.9261    999  0.026 *
Residuals 30 1.12892 0.037631                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(func ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(func))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.121064	0.361186	16.96203	0.001
Residual	30	1.982777	0.638814	NA	NA
Total	31	3.103841	1.000000	NA	NA

4.2.3 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "1_Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

antibiotics_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point == "2_Antibiotics") %>% 
  dplyr::select(sample) %>% pull()

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% #compute for each row (genome)
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
      !all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 

# A tibble: 9 × 2
# Rowwise: 
  phylum               sample_count
  <chr>                       <int>
1 p__Bacillota_A                 43
2 p__Bacteroidota                15
3 p__Bacillota                    8
4 p__Pseudomonadota               8
5 p__Cyanobacteriota              3
6 p__Verrucomicrobiota            2
7 p__Bacillota_B                  1
8 p__Bacillota_C                  1
9 p__Fusobacteriota               1

Antibiotics

structural_zeros %>% 
  filter(present=="antibiotics") %>% 
  select(1,5:10)

# A tibble: 4 × 7
# Rowwise: 
  genome                phylum            class                  order                 family                genus         species            
  <chr>                 <chr>             <chr>                  <chr>                 <chr>                 <chr>         <chr>              
1 AH1_2nd_7:bin_000003  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales   f__Enterobacteriaceae g__Proteus    s__Proteus cibarius
2 LI1_2nd_7:bin_000001  p__Bacillota_A    c__Clostridia          o__Clostridiales      f__Clostridiaceae     g__Sarcina    s__                
3 AH1_2nd_7:bin_000055  p__Bacillota      c__Bacilli             o__Mycoplasmatales    f__Mycoplasmoidaceae  g__Ureaplasma s__                
4 AH1_2nd_13:bin_000025 p__Bacillota_A    c__Clostridia          o__Christensenellales f__UBA3700            g__           s__                

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)

# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>

4.2.4 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output, 
     file = "data/ancombc_antib_accl.Rdata")

load("data/ancombc_antib_accl.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_time_point2_Antibiotics, p_time_point2_Antibiotics) %>%
  filter(p_time_point2_Antibiotics < 0.05) %>%
  dplyr::arrange(p_time_point2_Antibiotics) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_point2_Antibiotics)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point2_Antibiotics, y=forcats::fct_reorder(genome,lfc_time_point2_Antibiotics), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

Phyla of the significant MAGs

ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

            phylum sample_count
1     Bacteroidota           14
2      Bacillota_A            4
3        Bacillota            2
4 Campylobacterota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  select(phylum, genus, species)

             phylum           genus species
1  Campylobacterota  Helicobacter_J        
2       Bacillota_A                        
3      Bacteroidota     Bacteroides        
4      Bacteroidota     Odoribacter        
5      Bacteroidota     Bacteroides        
6         Bacillota   Coprobacillus        
7         Bacillota                        
8      Bacteroidota     Phocaeicola        
9      Bacteroidota     Bacteroides        
10     Bacteroidota     Phocaeicola        
11     Bacteroidota     Odoribacter        
12     Bacteroidota     Bacteroides        
13     Bacteroidota     Bacteroides        
14     Bacteroidota Parabacteroides        
15     Bacteroidota                        
16     Bacteroidota Parabacteroides        
17      Bacillota_A                        
18     Bacteroidota       Alistipes        
19      Bacillota_A                        
20      Bacillota_A                        
21     Bacteroidota     Odoribacter        

ancombc_rand_table_mag%>%
  filter(lfc_time_point2_Antibiotics<0)  %>% 
  count(genus, name = "sample_count") %>%
  arrange(desc(sample_count))

            genus sample_count
1                            6
2     Bacteroides            5
3     Odoribacter            3
4 Parabacteroides            2
5     Phocaeicola            2
6       Alistipes            1
7   Coprobacillus            1
8  Helicobacter_J            1

4.2.5 Differences in the functional capacity

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
  filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
    select(one_of(c("genome",sample_sub$sample))) %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)

4.2.5.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.331 0.0319
2 Antibiotics 0.244 0.134 

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94311, p-value = 0.09179

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Antibiotics)%>% 
  mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=treatment_colors) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

4.2.5.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment     MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.329 0.0319
2 Antibiotics 0.255 0.126 

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.95731, p-value = 0.2316

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))

Code_function	Acclimation	Antibiotics	Function	Difference	treatment
B06	0.68059020	0.47329940	Organic anion biosynthesis	0.20729080	Acclimation
B02	0.59930820	0.41576970	Amino acid biosynthesis	0.18353850	Acclimation
D02	0.38530780	0.20616790	Polysaccharide degradation	0.17913990	Acclimation
S03	0.27145162	0.09403129	Spore	0.17742033	Acclimation
B01	0.84159250	0.69078390	Nucleic acid biosynthesis	0.15080860	Acclimation
B07	0.44505530	0.30423320	Vitamin biosynthesis	0.14082210	Acclimation
B08	0.44618700	0.32094170	Aromatic compound biosynthesis	0.12524530	Acclimation
D09	0.25519710	0.13446900	Antibiotic degradation	0.12072810	Acclimation
B04	0.54481600	0.42941310	SCFA biosynthesis	0.11540290	Acclimation
D03	0.29173710	0.20687250	Sugar degradation	0.08486460	Acclimation
D05	0.28853330	0.22270070	Amino acid degradation	0.06583260	Acclimation
B10	0.05960097	0.03635236	Antibiotic biosynthesis	0.02324861	Acclimation

4.3 Does antibiotic administration remove the differences found in the two populations?

4.3.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(time_point == "2_Antibiotics" )

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(time_point == "2_Antibiotics" ) %>% 
  ggplot(aes(y = value, x = Population, color=Population, fill=Population)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
      facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

4.3.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(time_point == "2_Antibiotics") %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(time_point == "2_Antibiotics")

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$Population) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.015319 0.0153186 6.8764    999  0.021 *
Residuals 21 0.046782 0.0022277                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.356644	0.1527052	3.784762	0.001
Residual	21	7.527428	0.8472948	NA	NA
Total	22	8.884073	1.0000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$Population) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.030536 0.0305358 3.8593    999  0.061 .
Residuals 21 0.166158 0.0079123                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(neutral ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.785669	0.2085055	5.532084	0.001
Residual	21	6.778468	0.7914945	NA	NA
Total	22	8.564137	1.0000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$Population) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.012041 0.012041 0.9898    999  0.332
Residuals 21 0.255459 0.012165                     

adonis2(phylo ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.8963254	0.1888758	4.889993	0.001
Residual	21	3.8492558	0.8111242	NA	NA
Total	22	4.7455811	1.0000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$Population) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)
Groups     1 0.01804 0.018038 0.439    999  0.495
Residuals 21 0.86285 0.041088                    

adonis2(func ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.0184644	0.0098378	0.2086464	0.704
Residual	21	1.8584184	0.9901622	NA	NA
Total	22	1.8768828	1.0000000	NA	NA

4.4 Are the microbial communities similar in both donor samples?

samples_to_keep <- sample_metadata %>%
  filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2"))%>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2"))

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2")) %>% 
  ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000151 0.0001514 0.0181    999  0.905
Residuals 13 0.108675 0.0083596                     

adonis2(richness ~ time_point,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.06781753	0.02960432	0.3965972	0.979
Residual	13	2.22298047	0.97039568	NA	NA
Total	14	2.29079799	1.00000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00067 0.0006698 0.0523    999  0.801
Residuals 13 0.16645 0.0128035                     

adonis2(neutral ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.09261373	0.03872219	0.523666	0.804
Residual	13	2.29913452	0.96127781	NA	NA
Total	14	2.39174825	1.00000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00371 0.0037095 0.3076    999  0.635
Residuals 13 0.15676 0.0120588                     

adonis2(phylo ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.009763176	0.01920889	0.2546063	0.853
Residual	13	0.498500276	0.98079111	NA	NA
Total	14	0.508263452	1.00000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.013298 0.0132975 1.5346    999  0.217
Residuals 13 0.112644 0.0086649                     

adonis2(func ~ time_point,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	-0.003102642	-0.01647707	-0.2107297	0.812
Residual	13	0.191403206	1.01647707	NA	NA
Total	14	0.188300564	1.00000000	NA	NA

4.5 Does the donor sample maintain the microbial community found in acclimation?

4.5.1 Alpha diversity

sample_metadata <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
    TRUE ~ treatment
  ))
samples_to_keep <- sample_metadata %>%
  filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant"))

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

4.5.2 Beta diversity

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.008956 0.0089559 1.5113    999  0.272
Residuals 22 0.130368 0.0059258                     

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.2644898	0.06349573	1.491617	0.103
Residual	22	3.9009835	0.93650427	NA	NA
Total	23	4.1654732	1.00000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000794 0.0007940 0.0884    999  0.798
Residuals 22 0.197689 0.0089859                     

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3163831	0.07591112	1.807234	0.059
Residual	22	3.8514266	0.92408888	NA	NA
Total	23	4.1678097	1.00000000	NA	NA

neutral %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
 # scale_color_manual(values = treatment_colors,labels=c("Cold_wet" = "Cold wet", "Hot_dry" = "Hot dry")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.002011 0.0020114 0.2582    999  0.693
Residuals 22 0.171393 0.0077906                     

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.0409768	0.05601243	1.305392	0.286
Residual	22	0.6905894	0.94398757	NA	NA
Total	23	0.7315662	1.00000000	NA	NA

phylo %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     1 0.000522 0.0005217 0.052    999  0.828
Residuals 22 0.220603 0.0100274                    

adonis2(func ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.08891682	0.2161114	6.065212	0.046
Residual	22	0.32252296	0.7838886	NA	NA
Total	23	0.41143978	1.0000000	NA	NA

func %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
  scale_color_manual(values = treatment_colors,labels=c("Cold_wet" = "Cold wet", "Hot_dry" = "Hot dry")) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

4.6 Is the donor sample microbiota different to recipient microbiota?

4.6.1 Alpha diversity

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant"))

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
  stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

4.6.2 Beta diversity

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.12878 0.128776 13.013    999  0.001 ***
Residuals 20 0.19792 0.009896                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.76810	0.2746704	7.573672	0.001
Residual	20	4.66907	0.7253296	NA	NA
Total	21	6.43717	1.0000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.071939 0.071939 5.5113    999  0.038 *
Residuals 20 0.261057 0.013053                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	2.109391	0.3208221	9.447366	0.001
Residual	20	4.465565	0.6791779	NA	NA
Total	21	6.574956	1.0000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.04724 0.047241 2.6056    999  0.111
Residuals 20 0.36261 0.018131                     

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3764016	0.2391132	6.285119	0.001
Residual	20	1.1977547	0.7608868	NA	NA
Total	21	1.5741563	1.0000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000158 0.0001578 0.0131    999    0.9
Residuals 20 0.240836 0.0120418                     

adonis2(func ~ treatment,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.07897425	0.1809165	4.417536	0.084
Residual	20	0.35754884	0.8190835	NA	NA
Total	21	0.43652308	1.0000000	NA	NA

4.7 Does FMT change the microbial community over time?

4.7.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(type == "Treatment" & treatment %in% c("Post-FMT1","Post-FMT2") ) 

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(type=="Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2" )) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  scale_color_manual(values=treatment_colors)+
  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness

Model_rich <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Model_rich)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(4.3157)  ( log )
Formula: richness ~ treatment + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   171.5    174.8    -81.7    163.5       13 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.72426 -0.71683 -0.08993  0.38000  1.84050 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0        0       
Number of obs: 17, groups:  individual, 9

Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)          3.9416     0.1772  22.247   <2e-16 ***
treatmentPost-FMT2   0.4080     0.2420   1.686   0.0918 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
trtmnP-FMT2 -0.732
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Neutral

Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
      AIC      BIC   logLik
  129.347 132.1792 -60.6735

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.521142 11.47471

Fixed effects:  neutral ~ treatment 
                      Value Std.Error DF  t-value p-value
(Intercept)        24.81671  4.241888  8 5.850393  0.0004
treatmentPost-FMT2 15.06479  5.589803  7 2.695048  0.0309
 Correlation: 
                   (Intr)
treatmentPost-FMT2 -0.701

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.75054530 -0.55659466 -0.07097174  0.51610612  1.27053535 

Number of Observations: 17
Number of Groups: 9 

Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  56.51921 59.35141 -24.25961

Random effects:
 Formula: ~1 | individual
        (Intercept)  Residual
StdDev:   0.7091595 0.8466237

Fixed effects:  phylogenetic ~ treatment 
                      Value Std.Error DF   t-value p-value
(Intercept)        4.192894 0.3867555  8 10.841199  0.0000
treatmentPost-FMT2 0.939941 0.4163443  7  2.257606  0.0585
 Correlation: 
                   (Intr)
treatmentPost-FMT2 -0.582

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.1453262 -0.5912067 -0.1412005  0.4845880  2.3194143 

Number of Observations: 17
Number of Groups: 9 

Functional

Model_funct <- lme(fixed = functional ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_funct)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC   logLik
  -32.8173 -29.9851 20.40865

Random effects:
 Formula: ~1 | individual
        (Intercept)   Residual
StdDev:   0.0619285 0.03018711

Fixed effects:  functional ~ treatment 
                       Value  Std.Error DF  t-value p-value
(Intercept)        1.4883350 0.02345764  8 63.44777  0.0000
treatmentPost-FMT2 0.0323233 0.01501285  7  2.15304  0.0683
 Correlation: 
                   (Intr)
treatmentPost-FMT2 -0.352

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.08983687 -0.66983137 -0.03874894  0.52553040  1.62217086 

Number of Observations: 17
Number of Groups: 9 

4.7.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.017295 0.0172950 2.8289    999  0.111
Residuals 15 0.091706 0.0061137                     

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3571397	0.07959561	1.297184	0.00390625
Residual	15	4.1297869	0.92040439	NA	NA
Total	16	4.4869265	1.00000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.009753 0.0097526 0.7091    999   0.44
Residuals 15 0.206312 0.0137541                     

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3157079	0.08117526	1.325203	0.00390625
Residual	15	3.5735051	0.91882474	NA	NA
Total	16	3.8892130	1.00000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.010945 0.010944 0.6593    999  0.529
Residuals 15 0.248985 0.016599                     

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.06393539	0.09448624	1.565182	0.0703125
Residual	15	0.61272811	0.90551376	NA	NA
Total	16	0.67666350	1.00000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00162 0.0016197 0.1158    999  0.694
Residuals 15 0.20977 0.0139849                     

adonis2(func ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(func))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.02837271	0.07925317	1.291123	0.3203125
Residual	15	0.32962828	0.92074683	NA	NA
Total	16	0.35800099	1.00000000	NA	NA

4.8 Do FMT change the microbial community compared to antibiotics and acclimation?

4.8.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(type == "Treatment" & treatment %in% c("Antibiotics", "Acclimation","Post-FMT1") ) 

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
  filter(type=="Treatment" & treatment %in% c( "Antibiotics","Acclimation", "Post-FMT1")) %>% 
  ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
  geom_jitter(width = 0.2, show.legend = FALSE) +
  geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
#  scale_color_manual(values=treatment_colors)+
#  scale_fill_manual(values=treatment_colors) +
  facet_wrap(. ~ metric, scales = "free") +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

Richness: Problems to calculate

Model_rich <- glmer.nb(richness ~ treatment + (1|individual), data = alpha_div_meta)
#summary(Model_rich)
emmeans(Model_rich, pairwise ~ treatment)

Neutral

Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
#summary(Model_neutral)
emmeans(Model_neutral, pairwise ~ treatment)

$emmeans
 treatment   emmean   SE df lower.CL upper.CL
 Acclimation   17.4 3.46  8    9.410     25.4
 Antibiotics    8.5 3.91  8   -0.516     17.5
 Post-FMT1     24.8 3.67  8   16.296     33.2

Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                  estimate   SE df t.ratio p.value
 Acclimation - Antibiotics     8.90 4.82 13   1.847  0.1935
 Acclimation - (Post-FMT1)    -7.36 4.62 13  -1.591  0.2840
 Antibiotics - (Post-FMT1)   -16.25 4.92 13  -3.301  0.0148

Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 3 estimates 

Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
#summary(Model_phylo)
emmeans(Model_phylo, pairwise ~ treatment)

$emmeans
 treatment   emmean    SE df lower.CL upper.CL
 Acclimation   5.56 0.514  8     4.38     6.75
 Antibiotics   4.33 0.581  8     2.99     5.67
 Post-FMT1     4.17 0.544  8     2.92     5.43

Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                  estimate    SE df t.ratio p.value
 Acclimation - Antibiotics    1.229 0.727 13   1.690  0.2456
 Acclimation - (Post-FMT1)    1.387 0.698 13   1.986  0.1549
 Antibiotics - (Post-FMT1)    0.158 0.744 13   0.212  0.9755

Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 3 estimates 

Functional

Model_funct <- lme(fixed = functional ~ treatment, data = alpha_div_meta,
               random = ~ 1 | individual)
#summary(Model_funct)
emmeans(Model_funct, pairwise ~ treatment)

$emmeans
 treatment   emmean     SE df lower.CL upper.CL
 Acclimation   1.49 0.0437  8     1.39     1.60
 Antibiotics   1.47 0.0496  8     1.36     1.58
 Post-FMT1     1.49 0.0464  8     1.38     1.60

Degrees-of-freedom method: containment 
Confidence level used: 0.95 

$contrasts
 contrast                  estimate     SE df t.ratio p.value
 Acclimation - Antibiotics  0.02516 0.0661 13   0.380  0.9238
 Acclimation - (Post-FMT1)  0.00591 0.0638 13   0.093  0.9953
 Antibiotics - (Post-FMT1) -0.01924 0.0679 13  -0.283  0.9569

Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 3 estimates 

4.8.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "Post-FMT1")) %>% 
  select(sample) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "Post-FMT1"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.015456 0.0077280 1.0306    999    0.4
Residuals 21 0.157471 0.0074986                     

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.837434	0.201293	2.646248	0.001
Residual	21	7.290723	0.798707	NA	NA
Total	23	9.128157	1.000000	NA	NA

pairwise.adonis(richness, subset_meta$treatment, perm = 999)

                       pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 0.8277738 2.288600 0.1405032   0.001      0.003   *
2   Acclimation vs Post-FMT1  1 0.9551308 2.926054 0.1632291   0.001      0.003   *
3   Antibiotics vs Post-FMT1  1 0.9744543 2.741152 0.1741392   0.003      0.009   *

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.017834 0.0089172 0.5922    999   0.59
Residuals 21 0.316187 0.0150565                     

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	2.356335	0.2697712	3.879056	0.001
Residual	21	6.378232	0.7302288	NA	NA
Total	23	8.734566	1.0000000	NA	NA

pairwise.adonis(neutral, subset_meta$treatment, perm = 999)

                       pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 0.8812272 2.745592 0.1639591   0.001      0.003   *
2   Acclimation vs Post-FMT1  1 1.3127773 4.629826 0.2358567   0.001      0.003   *
3   Antibiotics vs Post-FMT1  1 1.3423457 4.351968 0.2508054   0.001      0.003   *

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.15470 0.077349 2.6532    999  0.098 .
Residuals 21 0.61221 0.029153                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.326050	0.3645754	6.024383	0.001
Residual	21	2.311195	0.6354246	NA	NA
Total	23	3.637246	1.0000000	NA	NA

pairwise.adonis(phylo, subset_meta$treatment, perm = 999)

                       pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 0.6891376 5.128491 0.2681074   0.002      0.006   *
2   Acclimation vs Post-FMT1  1 0.2531943 3.273842 0.1791546   0.036      0.108    
3   Antibiotics vs Post-FMT1  1 1.0996467 9.041595 0.4102060   0.003      0.009   *

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.23450 0.117249 4.3966    999  0.025 *
Residuals 21 0.56003 0.026668                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(func ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(func))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.636364	0.5726008	14.0672	0.002
Residual	21	1.221410	0.4273992	NA	NA
Total	23	2.857774	1.0000000	NA	NA

pairwise.adonis(func, subset_meta$treatment, perm = 999)

                       pairs Df  SumsOfSqs   F.Model        R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics  1 1.11238376 14.789477 0.5137112   0.001      0.003   *
2   Acclimation vs Post-FMT1  1 0.05911164  2.624236 0.1488993   0.137      0.411    
3   Antibiotics vs Post-FMT1  1 1.36464611 16.864488 0.5647004   0.001      0.003   *

4.9 Is the gut microbiota similar to the donor after FMT?

4.9.1 Beta diversity

sample_metadata <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
    TRUE ~ treatment
  ))

samples_to_keep <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>% 
  select(sample) %>% 
  pull()

subset_meta <- sample_metadata %>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.11342 0.113420 11.957    999  0.001 ***
Residuals 28 0.26559 0.009485                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.182046	0.154139	5.102365	0.001
Residual	28	6.486654	0.845861	NA	NA
Total	29	7.668700	1.000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.04487 0.044865 3.2331    999  0.079 .
Residuals 28 0.38855 0.013877                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(neutral ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.146508	0.1607364	5.362581	0.001
Residual	28	5.986340	0.8392636	NA	NA
Total	29	7.132848	1.0000000	NA	NA

neutral %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
  geom_point(size = 4) +
 # scale_color_manual(values = treatment_colors) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme(
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
    panel.background = element_blank(),
    axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 12),
    legend.position = "right", legend.box = "vertical"
    )

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     1 0.00008 0.0000787 0.005    999  0.937
Residuals 28 0.43669 0.0155959                    

adonis2(phylo ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.1296962	0.1018031	3.173567	0.002
Residual	28	1.1442941	0.8981969	NA	NA
Total	29	1.2739903	1.0000000	NA	NA

Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00120 0.001202 0.1035    999  0.744
Residuals 28 0.32533 0.011619                     

adonis2(func ~ treatment,
        data = subset_meta %>% arrange(match(sample,labels(func))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.01544812	0.02751097	0.7920986	0.46
Residual	28	0.54607759	0.97248903	NA	NA
Total	29	0.56152571	1.00000000	NA	NA

4.9.2 Structural zeros

FMT_samples <- sample_metadata %>%
  filter(type == "Treatment" & treatment == "FMT") %>% 
  dplyr::select(sample) %>%
  pull()

Transplant_samples <- sample_metadata %>%
  filter(type == "Treatment" & treatment =="Transplant") %>% 
  dplyr::select(sample) %>% pull()

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% #compute for each row (genome)
   mutate(all_zeros_FMT = all(c_across(all_of(FMT_samples)) == 0)) %>% 
   mutate(all_zeros_Transplant = all(c_across(all_of(Transplant_samples)) == 0)) %>% 
   mutate(average_FMT = mean(c_across(all_of(FMT_samples)), na.rm = TRUE)) %>% 
   mutate(average_Transplant = mean(c_across(all_of(Transplant_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_FMT == TRUE || all_zeros_Transplant==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_FMT & !all_zeros_Transplant ~ "Transplant",
      !all_zeros_FMT & all_zeros_Transplant ~ "FMT",
      !all_zeros_FMT & !all_zeros_Transplant ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "FMT", average_FMT, average_Transplant)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Structural zeros in each group

fmt <- structural_zeros %>% 
  filter(present=="FMT") %>% 
  count(phylum, name = "FMT") %>%
  arrange(desc(FMT)) 

fmt_f <- structural_zeros %>% 
  filter(present=="FMT") %>% 
  count(family, name = "FMT") %>%
  arrange(desc(FMT)) 

structural_zeros %>% 
  filter(present=="Transplant") %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant)) %>% 
  full_join(.,fmt, by="phylum" ) %>%
  mutate(across(everything(), ~ replace_na(., 0))) %>%
  as.data.frame() %>% 
  summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE)))

  Transplant FMT
1         52  55

Phyla to which the structural zeros belong in each group

structural_zeros %>% 
  filter(present=="Transplant") %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant)) %>% 
  full_join(.,fmt, by="phylum" ) %>%
  mutate(across(everything(), ~ replace_na(., 0))) %>% 
  tt()

phylum	Transplant	FMT
p__Bacillota_A	20	21
p__Bacillota	16	10
p__Pseudomonadota	6	11
p__Bacteroidota	3	6
p__Desulfobacterota	2	2
p__Bacillota_B	1	0
p__Chlamydiota	1	0
p__Cyanobacteriota	1	0
p__Fusobacteriota	1	0
p__Verrucomicrobiota	1	2
p__Actinomycetota	0	1
p__Bacillota_C	0	1
p__Campylobacterota	0	1

Families to which the structural zeros belong in each group

structural_zeros %>% 
  filter(present=="Transplant") %>% 
  count(family, name = "Transplant") %>%
  arrange(desc(Transplant)) %>% 
  full_join(.,fmt_f, by="family" ) %>%
  mutate(across(everything(), ~ replace_na(., 0))) %>% 
  tt()

family	Transplant	FMT
f__Lachnospiraceae	7	9
f__Erysipelotrichaceae	6	5
f__UBA660	6	0
f__Enterobacteriaceae	5	2
f__CAG-508	3	0
f__Ruminococcaceae	3	3
f__Anaerovoracaceae	2	0
f__Coprobacillaceae	2	0
f__Desulfovibrionaceae	2	2
f__Oscillospiraceae	2	1
f__Tannerellaceae	2	1
f__UBA1242	2	0
f__	1	3
f__Akkermansiaceae	1	0
f__CAG-239	1	2
f__Chlamydiaceae	1	0
f__Enterococcaceae	1	3
f__Fusobacteriaceae	1	0
f__Gastranaerophilaceae	1	0
f__Marinifilaceae	1	0
f__Mycoplasmoidaceae	1	1
f__Peptococcaceae	1	0
f__Anaerotignaceae	0	2
f__Bacteroidaceae	0	2
f__LL51	0	2
f__UBA3700	0	2
f__Acutalibacteraceae	0	1
f__Arcobacteraceae	0	1
f__CAG-274	0	1
f__CALVMC01	0	1
f__Devosiaceae	0	1
f__Mycobacteriaceae	0	1
f__Pumilibacteraceae	0	1
f__RUG11792	0	1
f__Rhizobiaceae	0	1
f__Rikenellaceae	0	1
f__Sphingobacteriaceae	0	1
f__Streptococcaceae	0	1
f__UBA1997	0	1
f__UBA3830	0	1
f__Weeksellaceae	0	1

4.9.2.1 Functional capacities of the structural zeros

#Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% structural_zeros$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

#Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

#Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
    TRUE ~ treatment
  ))%>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
  filter(genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",sample_sub$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
  column_to_rownames(., "genome") 
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=FMT-Transplant)%>% 
  mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=treatment_colors) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Elevation")

4.9.3 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
    TRUE ~ treatment
  ))%>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>% 
  column_to_rownames("sample") %>%
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "treatment",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output, file = "data/ancombc_FMT_transplant.Rdata")

load("data/ancombc_FMT_transplant.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output$res %>%
  dplyr::select(taxon, lfc_treatmentTransplant, p_treatmentTransplant) %>%
  filter(p_treatmentTransplant < 0.05) %>%
  dplyr::arrange(p_treatmentTransplant) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_treatmentTransplant)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

Differential abundance MAGs in each group

fmt_count <- ancombc_rand_table_mag %>%
  filter(lfc_treatmentTransplant<0)  %>% 
  count(phylum, name = "FMT") %>%
  arrange(desc(FMT)) 

ancombc_rand_table_mag %>%
  filter(lfc_treatmentTransplant>0)  %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant))  %>% 
  full_join(.,fmt_count, by="phylum")%>%
  mutate(across(where(is.numeric), ~ replace_na(., 0)))

             phylum Transplant FMT
1      Bacteroidota         13   5
2       Bacillota_A          4  18
3         Bacillota          3   1
4  Campylobacterota          1   1
5  Desulfobacterota          1   2
6     Spirochaetota          1   0
7    Pseudomonadota          0   5
8   Cyanobacteriota          0   2
9       Bacillota_B          0   1
10      Bacillota_C          0   1

ancombc_rand_table_mag %>%
  filter(lfc_treatmentTransplant>0)  %>% 
  count(phylum, name = "Transplant") %>%
  arrange(desc(Transplant))  %>% 
  full_join(.,fmt_count, by="phylum")%>%
  mutate(across(where(is.numeric), ~ replace_na(., 0))) %>% 
  as.data.frame() %>% 
  summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE)))

  Transplant FMT
1         23  36

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_treatmentTransplant, y=forcats::fct_reorder(genome,lfc_treatmentTransplant), fill=phylum)) + #forcats::fct_rev()
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
#  coord_flip()+
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

4.9.4 Differences in the functional capacities

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

sample_sub <- sample_metadata %>%
    mutate(treatment=case_when(
    treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
    TRUE ~ treatment
  ))%>%
  filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))

genome_counts_row <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
    select(one_of(c("genome",sample_sub$sample))) %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)

4.9.4.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment    MCI     sd
  <chr>      <dbl>  <dbl>
1 FMT        0.354 0.0255
2 Transplant 0.351 0.0426

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94866, p-value = 0.1557

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 128, p-value = 0.4826
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(treatment) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=FMT-Transplant)%>% 
  mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) 

means_gift <- element_gift_filt %>% 
  select(-treatment) %>% 
  pivot_longer(!sample, names_to = "Elements", values_to = "abundance") %>% 
  left_join(sample_metadata, by=join_by(sample==sample)) %>% 
  group_by(treatment, Elements) %>%
  summarise(mean=mean(abundance))

log_fold <- means_gift %>%
  group_by(Elements) %>%
  summarise(
    logfc_fmt_transplant = log2(mean[treatment == "FMT"] / mean[treatment == "Transplant"])
    ) %>% 
  left_join(., difference_table, by="Elements")

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
#  geom_point(size=4) + 
  scale_fill_manual(values=treatment_colors) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="Treatment")

log_fold %>%
  filter(Elements!="D0611") %>% 
  ggplot(aes(x=forcats::fct_reorder(Function,logfc_fmt_transplant), y=logfc_fmt_transplant, fill=group_color)) + 
  geom_col() +
  scale_fill_manual(values=treatment_colors) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Log_fold")+
  labs(fill="Treatment")

4.9.4.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(treatment) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  treatment    MCI     sd
  <chr>      <dbl>  <dbl>
1 FMT        0.349 0.0221
2 Transplant 0.354 0.0374

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.87186, p-value = 0.001843

wilcox.test(value ~ treatment, data=MCI)

    Wilcoxon rank sum exact test

data:  value by treatment
W = 121, p-value = 0.6801
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(12,13)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
    pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))